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HIGH  RESOLUTION  SIMULATION  OF  TURBULENT  FLOW 

IN  A  CHANNEL 


INTRODUCTION 

The  primary  purpose  of  the  current  work  is  to  perform  a  high  resolution,  statistically  stationary 
simulation  of  turbulent  flow  in  a  channel.  High  resolution  is  of  critical  importance  as  the  following 
review  of  the  experimental  facts  will  show.  It  is  well  known,  for  example,  that  the  turbulent  struc¬ 
tures  near  the  wall  are  highly  elongated  with  a  streamwise  wave  length  several  times  greater  than  the 
spanwise  wavelength.  The  characteristic  spanwise  wavelength  is  known  to  be  on  the  order  of  100 
viscous  lengths.  The  viscous  length,  /*,  is  defined  by  v/u* ,  where  v  is  the  kinematic  viscosity  and 
u*  is  the  friction  velocity  given  by  Vr^/p,  where  rw  is  the  shear  stress  at  the  wall  and  p  is  the  den¬ 
sity.  It  has  been  shown  by  numerous  experimental  studies  that  this  “streaky”  structure  is  clearly 
associated  with  a  process  called  bursting.  It  is  this  bursting  process  which  is  apparently  responsible 
for  the  production  of  new  turbulence  in  shear  flows.  For  this  reason,  it  is  essential  that  the  spanwise 
resolution  of  the  calculation  be  as  high  as  possible.  In  the  present  Cray-XMP  calculation  we  used  the 
highest  possible  resolution  consistent  with  our  code  and  with  current  machine  memory  limitations. 
The  simulation  represents  turbulence  in  a  channel  using  (128)  x  64  x  65  grid  points,  where  the 
spanwise  Ocj)  direction  was  represented  by  128  points,  the  streamwise  (x ,)  direction  by  64  points,  and 
the  wall-normal  ( x2 )  direction  by  65  points.*  The  channel  size  (based  on  channel  half-width)  is  2  in 
the  wall-normal  direction  and  5  in  both  the  streamwise  and  spanwise  directions.  This  gives  a  span- 
wise  resolution,  A xj/l*,  of  6.6,  which  should  be  more  than  adequate  to  resolve  the  near-wall  struc¬ 
tures. 
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Achieving  a  statistically  stationary  state  is  also  of  considerable  importance.  In  previous  direct 
simulations  of  turbulent  channel  flow  by  Orszag  &  Patera/ 11  the  flow  is  driven  by  a  constant  driving 
force  which  is  chosen  to  maintain  the  original  Poiseuille  flow.  The  introduction  of  highly  unstable 
disturbances  causes  transition  to  turbulence  so  that  the  wall  shear  stress  coefficient  increases  and  a 
quasistationary  state  is  achieved.  However,  the  mean  wall  shear  stress  must  eventually  return  to  its 
original  laminar  value  and  the  final  state  of  the  flow  may  or  may  not  be  laminar.  This  kind  of  simu¬ 
lation  is  unsatisfactory  since  it  produces  only  a  very  short  burst  of  quasi-stationary  turbulence.  To 
overcome  this  difficulty,  the  present  simulation  was  carried  out  in  two  stages.  First,  a  lower  resolu¬ 
tion  calculation  using  64  x  64  x  65  grid  points  was  started  from  laminar  initial  conditions  with  a 
driving  force  which  was  constantly  adjusted  to  match  the  wall  shear  stress.  This  resulted  in  a  reason¬ 
ably  stationary  flow  with  a  driving  force  of  about  6.06//?,  where  R  =  — .  This  initial  calculation 

V 

was  terminated  at  about  75  time  units,  where  time  is  made  non-dimensional  on  the  channel  half  width, 
h,  and  the  laminar  centerline  velocity,  U.  Once  this  quasistationary  state  was  achieved,  the  data  was 
interpolated  to  achieve  the  higher  spanwise  resolution  and  the  second  stage  of  the  calculation  was  ini¬ 
tiated.  At  this  point  the  driving  force  was  set  equal  to  6.0//? .  This  value  is  maintained  throughout 
the  rest  of  the  calculation. 

NUMERICAL  METHODS 

The  numerical  methods  used  here  have  been  developed  by  Orszag  and  others.  The  simulation 
involves  obtaining  approximate  solutions  to  the  following  problem: 

— -  =  u  X  3  —  V  n  +  -  V 2  u  +  / 
dt  R  J 

n  =  p  +  y  |« |2  0) 

v  ■  a  =  o 

•(The  notation  X  ,y ,  and  '  will  be  used  to  denote  streamwise,  spanwise.  and  wall-normal  coordinates  when  the  use  of 
subscripts  l.  3,  and  2  is  cumbersome.) 
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u  =  0  at  i,  =  ±1 


where  the  channel  walls  are  at  x2  =  ±1.  All  variables  in  equation  (1)  are  made  non-dimensional 
with  h,  U  and  p.  We  also  define  u,  the  velocity,  3,  the  vorticity,  p .  the  pressure,  and  /,  an  exter¬ 
nally  applied  force  field. 

The  problem  is  solved  by  a  pseudospectral  or  collocation  method  in  which  the  velocity  is 
expressed  by: 


M/2-1  V/2-1  J 

u(x,y  ,z ,/)  =  £  EL  smn,(t)e 

m»-A#/2  j  *0 


-  .  mx  t  ny 

2«  —  *  T 


Here,  X  and  Y  represent,  respectively,  the  streamwise  and  spanwise  lengths  of  the  domain.  The  prob¬ 
lem  is  discretized  by  choosing  the  collocation  points  as  follows: 


*,  =  i  -j 7  i  *  0...  M  -  I 

A# 


x-Jj j 


Zt  =  cos  k  -  0. .  .J  . 

J 

The  Fourier  coefficients,  a^j •  are  computed  efficiently  using  one  dimensional  fast  Fourier 
transforms.  Derivatives  then  can  be  calculated  in  Fourier  space  and  transformed  back  to  real  space 
without  introducing  phase  errors. 

The  time  stepping  is  split  into  three  intermediate  problems.  In  the  first  step,  the  non-linear 
term,  (u  x  3),  is  computed  directly  in  real  space  and  used  to  correct  the  current  velocity  field  using 
an  Adams-Bashforth  explicit  time  stepping  scheme  of  second  order.  The  incompressibility  condition 
is  then  imposed  by  inverting  the  Poisson  equation  for  the  pressure  and  satisfying  the  inviscid  boun¬ 
dary  condition  that  the  normal  velocity  on  both  walls  is  zero.  In  the  final  step,  the  viscous  correction 
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is  added  together  with  the  no-slip  boundary  condition  at  the  walls.  Further  details  on  the  numerical 
scheme  can  be  found  in  Orszag  &  Kells. (2) 

RESULTS 


A.  Global  Properties 

The  global  properties  of  the  flow  are  shown  in  Figures  1  and  2.  In  Figure  1  the  Reynolds 
number,  /?*,(/?*  =  u»h/v)  is  plotted  against  computational  time.  The  entire  computational  run  time 
(0  -  54.4)  is  shown.  It  is  important  to  note  that  in  a  physical  experiment,  the  flow  is  driven  by  a 
pressure  gradient.  In  the  simulation,  however,  periodic  boundary  conditions  are  used  in  the  horizon¬ 
tal  plane  so  that  a  differential  pressure  forcing  cannot  be  prescribed.  This  requires  the  use  of  the 
external  driving  force,  /,  in  the  equation  of  motion.  We  note,  of  course,  that  /  is  a  completely  arbi¬ 
trary  function  of  space  and  time.  If  we  define  a  control  surface  to  be  that  which  encloses  the  entire 
computational  domain,  a  simple  momentum  balance  in  the  streamwise  direction  results  in: 


iL,  <•■>*- 


-2  <r0> 
pU2 


2<f  i>  h 


In  this  expression,  <  >  represents  averaging  in  the  horizontal  plane,  t0  is  the  instantaneous  shear 
stress  at  the  wall,  and  f  x  is  the  streamwise  force  per  unit  volume.  The  simple  form  of  this  expression 
is  due  to  the  use  of  periodic  boundary  conditions  which  insure  no  net  mass  flow  through  the  control 
surface  and  no  net  force  arising  from  normal  stresses  (pressure).  It  follows  from  this  that  for  a  glo¬ 
bally  steady  flow  the  time  rate  of  change  in  momentum  within  the  control  volume  must  be  zero  so 
that  lim  <  r0  >  =  h  lim  <  f  x>  .  This  can  be  written  more  conveniently  as: 

I  —  m  t  —  a» 
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state  value  of  R*  =  v6 R  which  is  173.2  for  R  =  5000  (a  value  which  was  fixed  throughout  the  cal¬ 
culation).  We  can  see  readily  from  Figure  1  that  at  the  end  of  calculation  R*  is  within  about  0.5%  of 
the  theoretical  value  and  is  never  more  than  5%  lower  than  theory  throughout  the  entire  calculation. 
Thus,  for  all  practical  purposes,  the  turbulence  has  achieved  a  dynamic  equilibrium.  It  is  also  impor¬ 
tant  to  note  that  this  state  cannot  decay  to  a  laminar  or  quasi-laminar  state  as  in  the  calculation  of 
Orszag  because  of  the  maintenance  of  a  high  streamwise  driving  force.  There  is  no  reason  to  believe 
that  the  turbulence  cannot  be  maintained  indefinitely.  We  should  also  note,  however,  that  scaling 
arguments  show  that  the  flow  evolves  over  a  viscous  time  period  of  order  R  so  that  for  true  stationar- 
ity  we  must  calculate  for  about  5000  time  units.  This  is  clearly  well  beyond  the  currently  available 
resources  so  that  we  must  be  content  with  something  that  is  not  completely  stationary. 

In  Figure  2  we  plot  the  average  streamwise  velocity  at  the  center  of  the  channel,  Uct%  against 
computational  time.  The  centerline  velocity  has  been  made  nondimensional  with  the  initial  laminar 
centerline  velocity,  U.  It  appears  to  be  slowly  decreasing  during  the  course  of  the  run.  This  result  is 
not,  however,  inconsistent  with  the  fact  that  the  imposed  driving  force  is  always  larger  than  the  wall 
shear  stress  since  the  force  balance  (Eq.  3)  shows  that  the  mass  flux  (not  necessarily  Ucl)  must  be 
increasing  with  time.  It  is  also  interesting  to  note  that  the  R*  time  record  of  Figure  1  shows  high  fre¬ 
quency  fluctuations  which  are  not  present  in  the  centerline  velocity  data  of  Figure  2.  This  is  most 
likely  a  real  effect  due  to  the  much  shorter  time  scales  near  the  wall. 

We  also  examine  the  stationarity  of  the  flow  by  observing  changes  in  the  mean  streamwise  velo¬ 
city  profile  and  the  root  mean  square  amplitude  of  the  turbulence  intensities  as  the  calculation 
proceeds.  In  each  of  the  succeeding  figures,  we  plot  results  at  times  of  0,  12.825.  25.65.  38.48.  and 
51.3.  In  Figure  3  we  plot  the  mean  streamwise  velocity  scaled  on  u,  against  wall  normal  distance 
x",  defined  as  x2//*.  At  t  =  51.3,  three  regions  of  the  profile  are  clearly  identified:  (a)  the 
viscous  sublayer  (x2  s  2),  (b)  a  buffer  layer  (2  <  x2  <  30)  and.  (c)  a  logarithmic  layer 
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(30  sii*  s  120).  There  is  also  a  core  region  extending  out  to  x"  =  173.  The  significance  of  this 
is  that  very  early  in  the  calculation  the  logarithmic  layer  was  not  at  all  evident.  The  logarithmic 
region  at  t  =  51.3  can  be  approximated  by  the  so-called  law  of  the  wall, 
<«!>/«*  =  fc  In  x"  +  c ,  where  k  =  3.28  and  c  =  6.32.  Most  experimental  results  for  channel 
flow  give  k  between  2.5  and  3  and  c  varies  from  4  to  6.  The  calculated  values  are  somewhat  higher 
for  both  constants. 

In  Figures  4-7,  we  plot  the  rms  magnitudes  of  u1(  u2,  u 3,  and  utu2  against  x2 .  The  most 
interesting  feature  in  Figure  4  is  the  distinct  shift  in  the  x"  location  of  the  maximum  of  <u}  >  / u* 
towards  the  wall  as  the  calculation  proceeds.  The  location  of  the  peak  shifts  from  about  40  at  the 
start  of  the  run  to  about  23  at  the  end.  This  is  an  encouraging  trend  since  it  is  known  from  the  experi¬ 
ments  of  Kreplin  and  Eckelmann(3)  at  about  the  same  Reynolds  number  that  the  peak  occurs  at  about 
13.  Though  the  calculated  location  of  the  peak  is  still  far  from  experiment,  the  trend  suggests  that 
more  run  time  may  bring  the  peak  even  closer  to  the  wall.  The  calculated  amplitude  of  V <u  f>  /u  * 
is  somewhat  lower  than  the  value  of  2.85  found  in  (3).  In  the  Appendix  (Figures  A1  and  A2)  we 
show  both  the  mean  streamwise  velocity  and  streamwise  turbulence  intensity  profiles  in  greater  detail. 
The  peak  value  of  \<u2  >/u*  is  about  0.96  at  the  end  of  the  calculation  and  occurs  at  about 
x"  —  100,  as  shown  in  Figure  5.  This  is  very  close  to  the  value  of  1.0  obtained  experimentally  by 
Clark(4)  at  x"  =  70.  The  V  <«32>/«"  profile  peaks  at  about  1.31  at  the  end  of  calculation  and 
occurs  at  x*  =  50  as  shown  in  Figure  6.  Clark  measured  a  peak  of  about  1 .25  at  x  *  =  25.  The 
<U\U2>  /u*2  profile  peaks  at  about  0.64  at  the  end  of  the  calculation  and  occurs  at  x*  =  50. 

As  an  aid  in  examining,  at  least  qualitatively,  the  structure  of  the  turbulence,  we  plot  in  Figure 
8  (a-d),  contours  of  u,.  in  horizontal  planes  at  several  distances  from  the  wall.  On  each  figure  we 
plot  two  contour  levels,  (u !  -  <  a  1  > )  -  U\'  and  (u ,  -  <  u ,  > )  -  1/4  « , '  where  u  { '  is  the  rms 
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magnitude.  We  have  chosen  to  plot  these  contour  levels  since  they  correspond  to  the  low  speed 
regions  of  the  flow  normally  observed  in  experimental  flow  visualization  studies. 

In  the  sublayer  (x“  =  1  86)  we  note  that  the  structure  is  distinctly  elongated  in  the  streamwise 
direction.  In  the  buffer  layer,  (jc *  =  10.1),  a  similar  structure  is  apparent  but  the  average  size  of  the 


structures  has  grown.  It  cannot  be  determined  definitely  from  these  plots  if  a  distinct  spanwise 
wavelength  corresponding  to  the  well  known  streaky  turbulent  structure  actually  exists.  This  needs  to 
be  determined  by  more  detailed  spectral  analysis  in  the  future.  However,  if  one  can  imagine  a 
characteristic  wavelength  it  is  certainly  not  far  from  the  accepted  value  of  about  100,  as  mentioned 
earlier.  At  the  bottom  edge  of  the  logarithmic  layer  (x“  =  44.8)  the  structures  have  grown  notice¬ 
ably  larger,  and  well  within  the  logarithmic  layer  (jc*  =  99.0)  the  structures  have  grown  to  the  size 
of  the  flow  domain. 


B.  Conditional  Sampling  Results 


In  addition  to  the  global  results,  the  data  generated  by  the  numerical  simulation  has  been  exam¬ 
ined  by  conditional  sampling  techniques.  The  object  of  conditional  sampling  is  to  detect  regions  in 
the  flow  which  are  of  interest.  The  second  quadrant  Reynolds  stress,  where  the  fluctuating  streamwise 
velocity  is  negative  and  the  normal  velocity  is  positive  is,  of  considerable  interest  since  it  represents 
an  ejection  of  fluid  from  the  wall  region.  Also  of  interest  are  shear  layers  which  exist  close  to  the 
wail.  These  are  detected  by  the  Variable  Interval  Time  Averaging  (VITA)  technique.  These  two 
conditional  sampling  techniques  detect  turbulent  bursts,  which  are  believed  to  be  the  major  dynamical 
event  responsible  for  the  creation  of  turbulent  kinetic  energy  and  momentum  transport  from  the  wall. 


Two  sets  of  conditional  sampling  experiments  have  been  performed  using  the  data  set  described 
in  the  report.  The  first  set  of  experiments  is  a  direct  continuation  of  the  work  started  in  Leighton'5'. 
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Using  the  second  quadrant  Reynolds  stress  conditional  sampling  scheme,  turbulent  structures  were 
detected  in  the  simulated  turbulent  channel  flow.  The  detections  were  performed  at  10.2  and  16.8 
viscous  lengths  from  the  wall,  and  the  events  were  averaged  in  the  x-z  plane. 

Figure  9  shows  some  typical  results  of  the  second  quadrant  Reynolds  stress  conditional  sam¬ 
pling.  The  first  contour  plot  is  the  streamwise  fluctuating  coherent  velocity  at  the  detection  time. 
The  second  plot  is  the  fluctuating  normal  velocity.  Due  to  the  choice  of  the  coordinate  system  in  the 
computations,  fluid  flowing  away  from  the  upper  wall  has  a  negative  value.  In  the  final  figure,  the 
second  quadrant  Reynolds  stress  is  presented  and  is  also  positive  due  to  the  coordinate  system. 

Several  conclusions  can  be  drawn  from  these  averages: 

(1)  They  are  consistent  with  the  results  obtained  from  the  numerical  simulation  generated  on  a 
Cray  1  at  Los  Alamos  National  Laboratory,  and  described  in  (5). 

(2)  They  show  the  event  developing  two  structures.  The  region  of  second  quadrant  Reynolds 
stress,  representing  fluid  that  is  being  ejected  from  the  wall,  is  initially  farther  away  from  the  wall 
than  the  detection  site.  As  the  time  approaches  the  detection  time,  the  region  of  second  quadrant  Rey¬ 
nolds  stress  moves  toward  the  wall.  After  detection  the  region  continues  to  move  toward  the  wall.  A 
second  region  of  second  quadrant  Reynolds  stress  appears  slightly  downstream  and  above  the  old 
region.  This  later  region  has  not  been  detected  in  earlier  work. 

(3)  The  events  detected  closer  to  the  wail  are  detected  two  to  four  viscous  time  units  later  than 
events  farther  from  the  wall,  indicating  that  the  region  ejecting  the  fluid  is  moving  toward  the  wall. 

This  work  is  continuing  and  a  model  of  the  turbulent  ejection  incorporating  this  information  has  not 
yet  been  developed. 


SUMMARY  AND  CONCLUSIONS 


A  high  resolution  simulation  of  turbulence  in  a  channel  is  performed  using  x  64  x  65  grid 
points.  The  calculation  achieves  a  resolution  of  Ay//*  =  6.6  which  should  be  sufficient  to  ade¬ 
quately  resolve  the  organized  structures  near  the  wall.  In  addition,  a  dynamic  equilibrium  is  achieved 
in  which  the  shear  forces  at  the  wall  are  approximately  balanced  by  a  constant  driving  force.  The 
flow,  however,  is  not  completely  stationary  as  evidenced  by  a  slowly  dropping  centerline  velocity  At 
the  end  of  the  calculation,  the  mean  velocity  profile  exhibits  a  well  defined  logarithmic  layer 

In  general,  the  magnitudes  of  peak  root-mean-square  values  of  the  three  components  of  the  tur¬ 
bulent  velocity  fluctuations  and  the  velocity  correlation  <«,«;>  agree  well  with  experiment.  The 
peaks  occur  somewhat  farther  from  the  wall  than  experimental  measurements,  but  as  the  calculation 
proceeds,  these  locations  appear  to  be  moving  closer  to  the  wall.  This  behavior  suggests  that  more 
computational  time  is  needed  to  achieve  true  stationarity  and  results  which  are  closer  to  experiment. 
Contour  levels  of  the  streamwise  velocity  at  various  locations  above  the  wall  show  a  clearly  defined 
elongated  structure  near  the  wall.  The  structure  becomes  close  to  the  scale  of  the  horizontal  domain 
in  the  logarithmic  layer.  In  addition,  conditional  sampling  of  results  indicate  clearly  defined  coherent 
events. 
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i  —  Mean  streamwise  velocity  versus  distance  from  the  wall.  □,  t  =  0;  O,  r  =  12.825;  A 
25.65;  +,  t  —  38.48;  x  ,  t  =51.3.  Note:  These  symbols  apply  to  all  succeeding  figures. 
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Fig.  9  —  The  ensemble-averaged  fluctuating  velocity  and  the  second  quadrant  Reynolds  stress,  nor¬ 
malized  using  m*  for  a  turbulent  burst  detected  using  the  second  quadrant  Reynolds  stress  detection 
scheme.  The  (low  was  sampled  at  r*  =  16.5  and  K  =  4,  the  threshold  constant.  Note  that,  due  to 
the  coordinate  system  used,  a  negative  normal  velocity  implies  the  ejection  of  fluid  away  from  the 
wall. 


